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A method is described by which a function defined on a cubic grid (as 

from a finite difference solution of a partial differential equation) can be 

resolved into spherical harmonic components at some fixed radius. This 

Cm I has applications to the treatment of boundary conditions imposed at radii 

larger than the size of the grid, following Abrahams, Rezzola, Rupright et 

al. [1] . In the method described here, the interpolation of the grid data to 

^ ' the integration 2-sphere is combined in the same step as the integrations to 

^ ■ extract the spherical harmonic amplitudes, which become sums over grid 

(^ , points. Coordinates adapted to the integration sphere are not needed. 
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k> ; 1. INTRODUCTION 

5h ' Abrahams, Rezzola, Rupright, et al. [1, 2, 3] used spherical harmonics to inter- 

polate between a 3D grid on which sources of gravitational radiation were modeled 
and a distant boundary on which outgoing wave boundary conditions were imposed. 
While applying this method to a simplified model problem (scalar gravity [4] ) , we 
have developed an accurate and efficient method for extracting spherical harmonic 
amplitudes from the data on a rectangular grid. Spherical harmonic amplitudes 
are in principle defined as 



$i™(r, t)^(h rne, cfMr, 9, 0, t) <fn (i) 
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by integrating over a sphere of radius r. Here $ is the scalar field whose spherical 
harmonic decomposition 



^{r,e,(b,t)=J2'^im{r,t)Yr{e,(b) (2) 



is desired. 

Our method combines in one step the interpolation from the rectangular grid to 
the sphere and the integration over the sphere. It is motivated by the remark that 
a surface integral is the derivative of a volume integral 

<^Ur,t) = ^^j'jJYrHr',0,cj,,t)d^^dr' . (3) 

Volume integrals are relatively straightforward on a rectangular grid. A finite differ- 
ence approximation to (3) gives an integral over a shell of finite thickness covering 
the sphere r = constant. Our full implementation described below has some addi- 
tional features. 

2. CONTINUUM THEORY 

In order to find spherical harmonic amplitudes <!>;„ (i?, t) at a given fixed radius 
R, and also their radial derivatives dr^im{'r)\r=R, we first consider a complete 3D 
"Fourier" decomposition of $(a;,y, z,t) in a shell of thickness 2A around R. As a 
set of orthonormal functions on this shell we take the functions 

Y^i^{r,e,4>)^Rr,{r)Yr{e,4>) , (4) 

where the i?„ are obtained from Legendre polynomials, 



which have the property 

r-R+A 

R,,(r)Rk(vydr = 6^k (6) 

fl-A 

as a consequence of the orthogonality of the Legendre polynomials 

»i 

P„(a;)Pfc(a;) dx ^ 25nk/{2n + 1) . (7) 

Then the functions Ynim are a complete orthonormal set of functions on the shell 
S = {R — A < r < R + A} using the inner product 

{f\9)= f{x,y,z)g{x,y,z)dxdydz . (8) 

and its associated norm || / |p= (/|/). When a function ^{x, y, z) is smooth enough 
that the generalized Fourier series converges pointwise on the sphere r = R, one 
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has 



$(i?, e, 0, t) = ^ $„z„i?„(i?)r," 



(9) 
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SO that the spherical harmonic amphtudes for $ at r = i?, i.e., the coefficients of 



y,", are 



*;m(i?)-^$«(mi?«(i?) 



(10) 



Since Rn{R) ~ for odd n, this sum in numerical applications will typically have 
only one or two terms, e.g., n = and possibly n ~ 2. 

The coefficients ^nim are calculated from the volume integrals that are evaluated 
numerically on the rectangular grid: 



so 



^ImiR) 



^nim =11 Y„i^{r, 9, <?!>)$(r, 9, 0) dx dy dz 



(11) 



Y,Rn{R)Rn{r) 



Yi"^{9,(j))^{r,9,(f>)dxdydz . (12) 



3. DISCRETIZATION 

When the scalar function ^{x, y, z) is known only at points a; on a cubic grid G 
with spacing Ax = Ay ~ Az = k, the inner product (8) for our generalized Fourier 
analysis is relaced by 



{f\9) = E /»ff(^) 



Wx 



(13) 
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where each grid point must be given some weight Wx. If a cube of volume k^ 
centered on the grid point lies entirely outside the shell S, the weight will be zero. 
If it lies entrely inside S, the weight will be its volume k^. For other grid points we 
found it overly complicated to calculate, for the associated cube, the volume which 
overlaps the shell R~A<r<R + A. We simply use the result for a grid point 
lying on one of the rectangular axes (i.e., we imagined the cube rotated so its edges 
were parallel or perpendicular to the radial direction). Thus we took 







if 
if 



-i?| > A 



(A + ifc — |i? — r|)A:^ otherwise 



|r-i?| < A-ifc 



(14) 



where r is the Euclidean distance from the origin, and considered only the case 
A:<2A. 

With the inner product modified by this discretization, the function set (4) will 
no longer be orthonormal. One will find that 



{Ya\Yb) ~ Gab ~ Gba 



(15) 
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is no longer a unit matrix. Here A and B arc index groups [nlm) so that Y^ 
abbreviates the Ynim{r^9^(j)) of (4). The mean square error (as measured by this 
new norm 13) in the representation 

$(x,2/,z)«^$'4yA(a;,2/,z) (16) 

A 

is now minimized by the choice 

$^ = ^G^^(rB|$) (17) 

B 

where G is the matrix inverse to Gab- This holds for any finite subset of the 
Ynim one chooses to employ. Note that if one defines "adjoint" harmonics by 

Y^ = Y.^bG''^ (18) 

B 

then these coefficients in (16) can be computed as 

$^ = (^4]$) . (19) 



Similarly, by defining 



Rim = Y. Rn{R)Y'''"'{r, 9, cjy) (20) 



in which the adjoint Y appear, not the original Ya of equation (4), one calculates 
directly the spherical harmonic amplitudes of equation (10) by evaluating the sum 

$;m-(i?;m|*) • (21) 

It is important to notice that, in typical applications following the model of [1], 
the grid and the sphere r = i? at which a spherical harmonic decomposition is used 
do not change while the function $ or other fields are dynamically evolved or are 
relaxed toward a stationary state. Then a possible strategy is to store in memory 
the quantities Rim{x)wx at each grid point within the integration shell so that the 
computation of spherical harmonic coefficients becomes 

^im{R) = Y. ^'™(^) ^Mx) (22) 

xeG 

where only $ changes at each evolution or relaxation step and numerous other 
complications such as the Gab, the Rn{r), and the Wx are buried in a large table 
which does not change at each iteration of the p.d.e., nor even at each run of the 
program. 

4. EXAMPLE 

As a test to illustrate the application of these ideas, Mr. Keith Watt has run 
a test of the recovery of known spherical harmonic coefficients by these means. A 
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TABLE 1 
Spherical harmonic amplitudes: $ input, B recovered. 



I 


m 


•Pim 


t>lm 


7o error 








9.000000 


8.999315 


-0.0076 


1 


-1 


8.000000 


8.000664 


0.0083 


1 





7.000000 


6.999812 


-0.0027 


1 


1 


6.000000 


5.998907 


-0.0182 


2 


-2 


5.000000 


5.000006 


0.0001 


2 


-1 


4.000000 


4.000457 


0.0114 


2 





3.000000 


3.001447 


0.0482 


2 


1 


2.000000 


2.000002 


0.0001 


2 


2 


1.000000 


0.999975 


-0.0025 



cubic grid has each x,y, z coordinate varying at intervals of fc = 0.2 in the range 
— 1.3 < x,y, z < +1.3. A known function <I> is evaluated at these 14"^ — 2744 grid 
points. The sphere on which the resolution into spherical harmonics was desired 
was specified by i? = 1.0 and the thickness of the shell approximating this sphere 
was fixed by A = 0.15 in the notations of equation (14). Thus nonzero weight 
Wg was assigned to grid points g at radii in the range 0.75 < r < 1.25, which 
includes about 800 points. The approximately 1720 points at r > 1.25 and the 
approximately 220 points at r < 0.75 thus play no role in this anaylsis. 

The choice A — (3/4)fc was motivated by the fact that, in a one dimensional 
analogue of this analysis, this choice leads to Simpson's rule for the numerical 
integration and diagonalizes the 3x3 matrix Gab when one uses only n < 2. 

The test problem used 

'i> = J2<i>im{r/RyYr (23) 

with < / < 2 for a total of nine terms using real valued spherical harmonics. 
The input values $;,„ and the recovered values Bim are shown in the table. The 
accuracy was less (0.1% rather than the typical 0.01% shown in the table) when the 
radial dependence was changed to (i?/r)'+^; but the function used for the table may 
be more like the behaviors expected in applications where the sphere in question is 
outside the strong source region but within a wavelength of the center. 
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